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Mixed-state dynamics in one-dimensional quantum lattice systems: 
a time-dependent superoperator renormalization algorithm. 
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We present an algorithm to study mixed-state dynamics in one-dimensional quantum lattice 
systems. The algorithm can be used, e.g., to construct thermal states or to simulate real time 
evolution given by a generic master equation. Its two main ingredients are (i) a superoperator 
renormalization scheme to efficiently describe the state of the system and (ii) the time evolving block 
decimation technique to efficiently update the state during a time evolution. The computational 
cost of a simulation increases significantly with the amount of correlations between subsystems but 
it otherwise depends only linearly on the system size. We present simulations involving quantum 
spins and fermions in one spatial dimension. 
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The most interesting quantum phenomena involve 
strongly correlated many-body systems, but studying 
such systems — a central task in the areas of condensed 
matter physics, quantum field theory and, since recent 
years, also quantum information science [l|, |2| — has too 
often proven a formidable challenge. Indeed, in quan- 
tum many-body theory only a few exact solutions are 
available, while most analytical approximations remain 
uncontrolled. As a consequence, numerical calculations 
are of great importance. But even these suffer from a 
severe computational obstacle: an exponential growth of 
degrees of freedom with the system size that renders the 
direct simulation of most quantum systems prohibitively 
inefficient. 

And yet, ingenious methods such as quantum Monte 
Carlo techniques [2j can be used to approximately eval- 
uate, e.g., certain ground state properties in quantum 
lattice models. In one dimensional lattices, strikingly 
accurate results for quantities such as ground state ener- 
gies and two-point correlators can be obtained by using 
White's density matrix renormalization group (DMRG) 
|4|, a technique that has dominated most numerical re- 
search in the field since its invention more than a decade 
ago. Generalizations of the DMRG have also yielded ac- 
curate low energy spectra 5] or allowed for the simulation 
of real time evolution for small times 6] . 

Recently the time evolving block decimation (TEBD) 
algorithm (jj has been proposed to simulate real time 
evolution in one-dimensional quantum lattice systems. 
This technique can be easily adapted to standard DMRG 
implementations |3, |9j and seems to be very efficient gj, 
9, 10]. As in DMRG, a decisive factor in the performance 
of the TEBD method is that not a lot of entanglement is 
present in the system, a condition that is ordinarily met 
in one-dimensional lattices at low energies |7|. 

In this paper we extend the TEBD algorithm to han- 
dle mixed states. We describe how to efficiently sim- 
ulate, in one-dimensional quantum lattice systems, real 
time Markovian dynamics as given by a (possibly time- 
dependent) master equation made of arbitrary nearest 



neighbor couplings. By considering evolution in imagi- 
nary time, the present extension can also be used to con- 
struct thermal states for any given temperature. Thus, 
we show how to numerically explore non-equilibrium 
many-body dynamics under realistic conditions, includ- 
ing the effects of finite-temperature and decoherence. 

A key observation for the success of the algorithm is 
that in one spatial dimension many states of interest, in- 
cluding thermal states and local perturbations thereof, 
contain only a restricted amount of correlations between 
subsystems, in a sense to be further specified. This fact, 
parallel to the restricted amount of entanglement ob- 
served in the pure-state case, allows us to introduce an 
efficient decomposition for the state of the system, re- 
ferred to as matrix product decomposition (MPD). The 
MPD is nothing but a mixed-state version of a matrix 
product state |ll|. and as such, we can use the TEBD 
to update it during a time evolution. It also follows that 
our scheme can again be fully incorporated into standard 
DMRG implementations without much programming ef- 
fort 0,0. 

We consider a generic one dimensional quantum lattice 
made of n sites, labeled by index /, / £ {1, • • • ,n}, each 
one described by a local Hubert space V) 1 ' = Cd of finite 
dimension d. We assume the evolution of the n sites, in 
a global state p, is given by a master equation |2| 

P = C\p] (1) 

where H and L^ are the Hamiltonian and Lindblad op- 
erators, and where we require that the (possibly time- 
dependent) Lindbladian superoperator C further decom- 
pose into terms involving at most two contiguous sites, 



c\p] = J2^+M- 



(2) 



Reduced superoperator s. — A pure state evolution is de- 
scribed by a vector |$) in the n-fold tensor product of 



Cd- Let us divide the n sites into two blocks, denoted L 
(left) and R (right). Then DMRG and TEBD consider 
reduced density matrices, e.g., that of block L, 



l*>e 



,M = 



tr«(|*X*l) €L(HW), (3) 



where L(H) denotes the set of linear mappings on 
H or, equivalently, the complex vector space of 
dim(H)xdim(H) matrices. Here we are concerned with 
the evolution of a mixed state, which requires more nota- 
tion. For each site I, let K^ = L(HW) = C d 2 denote the 
vector space of d x d complex matrices. We switch into 
representing a density matrix a G L(H) as a "superket" 
\a\ £K = L(H), while a superoperator Q_| L(L(H)) is 
regarded as a linear mapping Q G L(K) J12J . 



|$) G H 
cr G L(H) 
Q G L(L(H)) 



H e IK (4) 

Qt G L(K), 



where 1$), = ||$)($|) r For d x d matrices A and B, the 
scalar product B (|) between super kets |A). and |B) , and 
the action of Q on \A) V are defined through 



,(A|B), = itr(AtB) 



QaL4> =|QL4]>, 



(5) 



Also, if Q is a superoperator on a bipartite space rf L l <g> 
H^l and {jM^} is an orthonormal basis in K^ = 
L(H^l), we define the partial trace of Qt over block R as 




FIG. 1: Quantum Ising chain with transverse magnetic field, 
Eq. 1171 . at finite temperature, with local dimension d = 2, 
n = 100 sites and effective xt = 80. At zero temperature this 
model corresponds to a quantum critical point. The spectrum 
{A«^} of the reduced superoperator Qt^ L ' for the left n/2 = 50 
spins is plotted as a function of /3 g [0, 10/ J] (only the largest 
52 Aj a 's are shown). For any inverse temperature j3, a fast 
decay of the eigenvalues Aj^ in a ensures that the state can 
be accurately approximated by a MPD with small effective 
X»- 



tr.^QO^EW^U 



(6) 



Finally, let p G L(C d ®") be the state of the n-site lattice 
and |p) B its superket. We define the reduced superoperator 
for a block of sites, say for block L, as (see example |l3j) 

He(Q 2 ) 0n — a m =WIP>,(/>l) eL(K [L] ), (7) 

in analogy with Q, and rewrite equation Q) as 

\p\=C*\p\, (8) 

which parallels the Schrodinger equation |\I>) = —iH\fy). 

Renormalization of reduced superoperator s. — Given 

blocks L and R, the Schmidt decomposition of \p\ reads 

\p\ = Y, ^\Mi% ® \M l a R \, X la > X ta+1 > 0, (9) 

where the Schmidt super kets {\Ma ' )„} fulfill 

Q^\M [ a L \ = (A« Q ) 2 |AfI L] ) 8 , a m = tofldp^l), (10) 
Q,l*]|i\A = (A, a ) 2 |M^ ] ) t , Q„M = tHid^pD.Cll) 



The rank X« of the reduced superoperators Qt ' L ' and Q ^ R ' 
measures the amount of correlations between blocks L 
and R 14] . In principle its value is only bounded above 
by the dimensions of K\ L ^ and K\ R \ which grow expo- 
nentially in the number of sites. However, as the exam- 
ples below illustrate, many situations of interest involving 
one-dimensional mixed-state dynamics are only slightly 
correlated, in that the coefficients {A« Q } decay very fast 
with a. That is, a good approximation to \p) } can be 
obtained by truncating JUJl so that only a relatively small 
number of terms are considered. 

Thus, whereas an efficient description of |$) in J5J is 
achieved both in DMRG and TEBD by conveniently dec- 
imating, say, the block space H^ supporting the reduced 
density matrix p<- L ' , our efficient description of p is based 
instead on decimating the block space K^l = LHHI^l) 
supporting the reduced superoperator Q^ in (0) J15lll6j . 

Matrix Product Decomposition and TEBD. — We re- 
gard the n-site p as a vector \p\ in the ro-fold tensor 
product of C<p , while the master equation JSJ is formally 
identical to the Schrodinger equation. Simulating mixed- 
state dynamics can therefore be achieved by conveniently 
adapting the pure-state techniques of [jj. More specifi- 
cally, given an orthonormal basis {\ii\ } of K^ for each 



site I (I = 1, • • ■ , n), we expand \p\ as [l7j 



H = H-"H C ^- 



l*i>. ® 



(12) 



ii=0 



We choose |0z)„ = \I/d\ to be proportional to the identity 
(that is, as a mapping on E[W), so that physical normal- 
ization of p, tr(p) = 1, corresponds to cq.-.q = 1. We use 
a MPD for the coefficients c^...^, 

r . . V^ p[l]n\Jl]p[2]J2 \[2]-p[3]« 3 ...r[n]i„/ 1 o\ 



which can be built through a succession of Schmidt de- 
compositions of \p) f (see for details). Finally, we can 
use the TEBD method to update the tensors jT^i and 
{A«W} during an evolution of the form CQl-JU 18]. 

Example 1: Thermal state. — Given a nearest neighbor 
Hamiltonian H and an inverse temperature (3 = 1/kT, a 
mixed state of interest is the thermal state 



-0H 



Pfi 



Z(f3) Z(p) 



E« 



-0E, 



E S ){E S 



(14) 



where Z((3) = tr(e~P H ) is the partition function. One 
can numerically simulate pp by attempting to compute 
all relevant energy eigenstates \E S ) and averaging them 
with weights e~ /3Ea /Z((3). A very simple and efficient 
alternative is to build a MPD for \pp), by simulating 
an imaginary time evolution from the completely mixed 
state, 



-pH 



V=exp(-/3T)|/L 



(15) 



where super ket |/Y and superoperator T« correspond to 



in - \h\ ® 



\In\, T[A] = ~(HA- 



AH). (16) 



Indeed, exp(— [3T) can be Trotter expanded into transfor- 
mations involving only two adjacent sites, and the MPD 
can be therefore updated using the TEBD [lj|. Notice 
that a single run of the simulation builds the thermal 
state ppi for any intermediate value of (5' £ [O,0\. Fig. 
(|T|) corresponds to thermal states for a quantum Ising 
model with transverse magnetic field, 



* = £*? 



'i+i 



i-i 



i=i 



(17) 



Example 2: Time- dependent master equation. — We 
consider a lattice of n = 100 sites loaded with n/2 
fcrmions that evolve according to a Lindbladian 



C[p] = -i[H,p] +j^2(mpni 



;P n i 



nip), 



(18) 
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FIG. 2: Fermionic lattice of Eq. 118H at finite temperature 
j3 — 1/J, dephasing 7 = 0.4J, and with d = 2, n — 100. 
The particle current, see Eq. 1201 . is due to a time-dependent 
applied bias p(t) with turn-on time to = 2/J and rise time 
t s = 0.1/J. Simulations with an effective \t — 30, 40, 50 show 
rapid convergence to the exact solution. This convergence 
will be addressed in more detail in [l4||. 




FIG. 3: Same system as in Fig. ©• The spectrum {Aj Q } of 
the reduced superoperator Qj for the left n/2 = 50 sites is 
plotted as a function of time. The number of relevant eigen- 
values Att^ (say above 10 -6 ) increases as the applied bias is 
turned on, but remains small throughout the evolution and it 
even decreases for long times. 
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FIG. 4: Fermionic lattice of Eq. (I18H at finite temperature 
/3 = 1/J, dephasing 7 = 0.4J, and with d — 2, n — 100 and 
no applied bias, fi(t) = 0. Unequal time, two-point correlator 
t2"Tl for Z = 2, 4, 6, 8, 10 and t G [0, 10/ J]. The results corre- 
sponding to an effective \t — 40 and 50 practically overlap at 
all times, as the inset shows. 



ni = ojoj) where the last term accounts for phase damp- 
ing and the Hamiltonian part corresponds to hoping be- 
tween adjacent sites and a time-dependent on-site energy, 



n— 1 n/2 n 

H = -jJ2(a f ia l+1 +h.c.)-n{t){J2ni-J2 n i)> ( 19 ) 



z=i 



(=1 l=n/2+l 



where p(t) = po[e ~ + 1] _1 introduces a bias /Lto be- 
tween the left and right halves of the lattice at t = to- 
Fig. J5J shows the particle current 



-21m (oJo(t)o 5 i(t)> 



(20) 



from the right half of the lattice to the left half as a 
result of switching on the bias. The numerical results 
have been obtained by mapping the fermionic lattice into 
a spin lattice through a Jordan- Wigner transformation 
and by simulating the resulting spin lattice model, which 
still contains only nearest neighbor couplings, using the 
TEBD algorithm. The exact solution can be efficiently 
obtained by numerically integrating the time evolution 
of two point-correlators. Comparison with the exact 
solution shows that the simulations for small effective 
Xtt = 30, 40, 50 are remarkably accurate and converge for 
increasing value of Xt- The fast decay of the spectrum in 
Fig. J3J justifies this convergence. 

Example 3: Unequal-time correlators. — For the above 
fermion system with no bias, p,(t) = 0, and finite tem- 



perature, we finally consider the expectation value 

(aJ(i)ai(O)) = tr (a\£ t [a lP \) , (21) 



where St is the time evolution operator resulting from 
the master equation. The simulation (see Fig. J3J|) is 
achieved as follows: (i) the initial state of the system, 
a thermal state with (3 = 1/J, is obtained by evolution 
in imaginary time as explained in Example 1; (ii) the 
annihilation operator a\ is applied to the initial state p 
to obtain a\p; (Hi) a\p is evolved in time according to £t, 
(iv) the creation operator a \ is applied on £ t [a\p}] and (v) 
the trace of the resulting operator a\£t\aip\ is computed. 
Each of these steps can be performed efficiently by using a 
MPD and the update techniques of |jfl . In this particular 
case the Lindbladian C is time-independent and we can 
integrate Eq. JSJ), so that step (Hi) becomes 



\£t[aip}\ =£u\aip\ =exp(£.*)|aip), 



(22) 



Because of property J5J, exp(C t t) can be Trotter- 
expanded into small transformations involving only two 
adjacent sites, and therefore it can be implemented using 
the TEBD. 

We have presented an extension of the TEBD algo- 
rithm to mixed states. With specific examples involv- 
ing spins and non-interacting fermions, we have shown 
how to (i) construct thermal states; (ii) evolve a state in 
time according to a time-dependent master equation; and 
(Hi) compute unequal time correlation functions. The al- 
gorithm can be used for generic one-dimensional lattice 
systems, including interacting fermions and bosons as ex- 
plored in |l4| . 

MZ acknowledges support from an NSF Graduate Fel- 
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See also F. Verstreate et al., cond-mat/0406426 

Note: The extension of the TEBD method presented 
here was outlined in [G. Vidal, quant-ph/0301063 v2 
(2003)] but excluded from the first reference in jjj by 
request of a referee. 
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ijkl 
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rank(p [L1 ) = dim(H [L1 



\p\~\I m \®\I™\ 
rank(Q [L1 ) = 1. 
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\p\= J2 ^i^kiJj ® Wh\ G C 4 ®' 
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Qt = tm\p\( P \ (25) 

= ^2i< J i2\p\lp\° i A ( 26 ) 
»2 

= J2(J2 c i^2\^i\)(^2^i'i2»(^i'\) ( 27 ) 

*2 «1 «l' 

= ^2(^2ciivtCii'ia)Wii) t (' 7 ii'\ ( 28 ) 

= Z>;k4>.-i- ( 29 ) 



will have Xt — X coefficients given by Ajr 
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-0H 
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